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Abstract. 

Fluctuation relations for the entropy production in non equilibrium stationary 
states of Ising models are investigated by Monte Carlo simulations. Systems in contact 
with heat baths at two different temperatures or subject to external driving will be 
studied. In the first case, by considering different kinetic rules and couplings with the 
baths, the behavior of the probability distributions of the heat exchanged in a time 
t with the thermostats, both in the disordered and in the low temperature phase, 
are discussed. The fluctuation relation is always verified in the large r limit and 
deviations from linear response theory are observed. Finite-r corrections are shown to 
obey a scaling behavior. In the other case the system is in contact with a single heat 
bath but work is done by shearing it. Also for this system the statistics collected for 
the mechanical work shows the validity of the fluctuation relation and preasymptotic 
corrections behave analogously to the case with two baths. 
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1. Introduction 

In equilibrium statistical mechanics the knowledge of general expressions for the 
probabilities of microscopic configurations is the cornerstone of a successful theory 
describing the macroscopic behavior of a large variety of physical systems. Similar 
expressions, however, are not known for systems in non-equilibrium steady states 
(NESS), despite their widespread occurrence in nature and their practical interest. 

NESS are usually realized by driving a system, either mechanically, as in the case 
of sheared or stirred fluids, or thermodynamically, due for instance to couplings to 
reservoirs at different temperatures. These states are characterized by a finite rate 
of entropy production, and the recent proposal [H [2j [3] of a relation governing the 
fluctuations of this quantity constitutes an important result of general validity. The 
relation connects the probability P(S(r)) of producing an entropy S(r) in a time interval 
r, with the probability of the opposite quantity, according to 

Eq. ([1]), also known as Gallavotti-Cohen relation, holds in the large r limit, specifically 
with r larger than all relaxation times of the system. It was proved as a theorem 
for a specific class of dynamical systems in [3] and then established for stochastic 
kinetics in (U El E]. Expressions related to Eq. ([1]) have been established in Refs. 
El El EH EEl E21 CE3] • Recent reviews are given in [H] . 

Fluctuation relations (FRs) are expected to be relevant in mesoscopic systems, 
particularly in nano- and biological sciences |15j . because at these scales typical thermal 
fluctuations may be sufficiently large to be comparable with the magnitude of the 
external driving. FRs have nowadays been tested in some experiments [T6J, [171 HB] • 

In this work we study the FRs in simple standard statistical models, where their 
validity can be ascertained, and some of the mechanisms of their occurrence can be 
investigated. We consider the Ising model as a paradigmatic example of interacting 
system interested by a phase transition. The model is maintained in NESS either by 
coupling it to two thermal baths at different temperatures Ti,T 2 , or by a mechanical 
forcing. In the former case, introducing A/3* 1 ) = — A{3& = (^r — ^-J , the relation (PQ) 
can be specified as 

where Q^ n \r) is the heat exchanged with the heat bath at temperature T n [n = 1, 2) in 
a time r. With mechanical driving Eq. ([ID can be cast as 

P(W(t)) W(r) 
%(-W(r)) " T ' [6) 

where W(r) is the work done on the system in the time interval r and T is the 
temperature of the thermostat. 
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In the case of thermodynamic driving, the relation §Z§ has been previously studied 
for different systems. It was shown to hold [19] for a chain of oscillators coupled at the 
extremities to two thermostats, and it was studied [20] for a Brownian particle in contact 
with two reservoirs. A fluctuation relation has been proved in [21] for the heat exchanged 
between two systems initially prepared in equilibrium at different temperatures and later 
put in contact. The case of the Ising model coupled to different reservoirs, similar to 
that considered in this paper, has been studied analytically in [22], in a mean-field 
approximation. Fluctuation properties of work due to a magnetic field in transient 
regimes of Ising model have been analyzed in [23] . We are not aware of studies of the 
relation (J3l) in stationary states of mechanically driven Ising models. 

In this Article we show that the FRs hold in the large-r limit for the nearest 
neighbour Ising model and both kinds of NESS analyzed. In the case of NESS induced 
by the presence of two thermostats, we also consider different couplings with the baths 
and kinetic rules, some of which have been previously reported in |24j, in order to address 
the generality of the results. Our numerical data allow to appreciate and characterize 
the finite time corrections to the asymptotic result. These were shown to be of order 
1/r in 0,1251 EH]- For systems described by a Langevin equation, in cases corresponding 
to the experimental setup consisting of a resistor and a capacitor in parallel, finite time 
corrections also behave as 1/r when work fluctuations are considered [TTl ITT] , while 
faster decays have been predicted for other topologies of circuits [TTJ. Our data show 
that the leading term of such corrections decays as 1/r. We also propose an expression 
which well describes the corrections in an extended range of values of r, incorporating 
the sub-leading behavior. This expression implies a scaling behavior which takes into 
account the geometry of the system and the nature of the coupling with the heat baths. 

The occurrence of a phase transition in the Ising model allows us to discuss the 
interplay between the breaking of ergodicity and the validity of the fluctuation relation. 
A finite size system in the low temperature NESS remains trapped into broken symmetry 
states for a time r erg that diverges in the thermodynamic limit, much like in equilibrium. 
By varying the system size and the baths temperatures, we are able to investigate the 
large-r limit both in the regime r 3> r erg and r <C T erg . The latter is particularly 
interesting because in this case, since FRs are expected to hold for r much larger than 
the characteristic relaxation times of the systems, they are not necessarily obeyed in 
this condition. Interestingly, instead, we find that, the FR (jSJ) holds true in any case. 

Close to equilibrium the FR implies the Green-Kubo relation (GKR) |27j . 
Therefore, as discussed in [26] a stringent test of the FR, which cannot be reduced 
to linear response theory, can be only achieved when the drive is large enough to bring 
the system far from equilibrium, spoiling the GKR and/or determining non-Gaussian 
V(Q^). An analysis of the data in the case of thermodynamic driving will allow us to 
provide a strict test of the validity of the FR also in the far from equilibrium regime. 

The paper is organized as follows. In the next section we introduce the model and 
discuss two different implementations of the coupling to the heat baths. Then the results 
of our simulations will be presented and interpreted in terms of scaling expressions for 
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finite time corrections. Relations with linear response theory are discussed in Sect. 2.3. 
In Section 3 the model with mechanical driving is considered. Section 4 includes our 
conclusions and a discussion of the perspectives of this work. 

2. Systems in contact with heat baths at two temperatures 

2.1. The models 

We consider a two-dimensional Ising model defined by the Hamiltonian H{a} = 
~JJ2(ij)< T i cr j> where = ±1 is a spin variable on a site i of a rectangular lattice 
with N = M x L sites, {a} is the configuration of all the spins and the sum is over all 
pairs (ij) of nearest neighbors. 

In the case of statical coupling with the heat baths the system is divided into 
two halves. The left part (the first M/2 vertical lines) interacts with the heat bath 
at temperature T\ while the right part is in contact with the reservoir at T 2 > T\. 
We have implemented both open or periodic boundary conditions. We used Monte 
Carlo spin-exchange (Kawasaki) dynamics, corresponding to systems with conserved 
magnetization. The case with single spin dynamics was considered previously [24] . In 
the present case with two heat baths, we have implemented the Kawasaki rule as follows: 
picking at random a couple of nearest neighbor spins 07, a m we attempt their exchange 
according to standard Metropolis transition rates 



where {a} and {a'} are the configurations before and after the move and 
AE({a}, {a'}) = H{a'} - H{a}. The temperature T to be entered in A({a'}, {a}) 
is chosen as follows: if both the spins considered are in contact with the same 
bath, T is the temperature of that thermostat. In the case in which, say, o m 
is coupled to the temperature T x and o~\ to T 2 , we compute AE({a}, {cr'})/T by 
splitting the two contributions from the different reservoirs, namely AE({a}, {a'})/T = 



-(J/T 2 )a' l J2{ i ) l °'i-(J/ T i)v'mY,{i) m ^ + (^/ T 2)^E(i) ; Vi + (J/ T i)vmY,(i) m °h where (i), 



((i) m ) are nearest neighbors of 07 (o~ m ). 

In the second implementation, using single-spin dynamics, each spin <7j, at a given 
time t, is put dynamically in contact with one or the other reservoir depending on the 
(time dependent) value of hi = (1/2) | ^2(j\.crj\, where the sum runs over the nearest 
neighbor spins o~j of 07 . Notice that hi is one half of the (absolute value) of the local 
field. In two dimensions, with periodic boundary conditions, the possible values of hi 
are hi — 0, 1, 2. At each time, spins with hi = 2 are connected to the bath at T = T\ 
and those with hi = 1 with the reservoir at T = T 2 . Namely, when a particular spin o~j 
is updated, the temperature T\ or T 2 is entered into the transition rate according to the 
value of hi. Loose spins with hi = can flip back and forth regardless of temperature 
because these moves do not change the energy of the system. Then, as in the usual Ising 
model, they are associated to a temperature independent transition rate. Notice that 
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hi = 2 correspond to spins whose surrounding neighbors are aligned, a situation which 
is typically found in the bulk of ordered domains, while hi = 1 corresponds to interfacial 
spins. This model was introduced in [28] and studied in [29]. It is characterized by a line 
of critical points in the plane T 1? T 2 , separating a ferromagnetic from a paramagnetic 
phase analogously to the equilibrium Ising model. Metropolis transition rates have been 
considered also in this case. 

Considering the stationary state, a generic evolution of the system is given by the 
sequence of configurations {o~(t)} = {o~i(t), . . . ,o"Ar(t)} where o~i(t) is the value of the 
spin variable at time t. Denoting with t£ the times (measured here as the number of 
elementary Montecarlo updates) at which an elementary move is attempted by coupling 
the system to the n-th reservoir, the heat released by the bath in a time window [s, s + r] 
is defined as 

QW(r) = W^)) " # H4 n) - I))]- (5) 

{t[ n) }c[s,s+r] 

In the case of dynamic coupling to the thermostats, recalling the discussion above, 
(r) and (r) will be also referred to as bulk and interface exchanged heats. The 
properties of Q^ n '(r) will be computed by collecting the statistics over different sub- 
trajectories obtained by dividing a long history of length tp into many {tp/r) time- 
windows of length r, starting from different s. 

Notice that all the dynamical rules considered insofar obey the generalized detailed 
balance condition [30] 

e S(1)/Tl+2(2)/T M({a'}, {a}) = A({a}, {a'}) (6) 

where QP~\ are the heats exchanged with the reservoirs during an elementary 
transition. 

2.2. Results for T U T 2 > T c 

We begin our analysis with the study of the relation (T5]) in the case with both 
temperatures above the critical value T c ~ 2.269 of the equilibrium Ising model. In 
the following we will measure times in montecarlo steps (MCS) (1 MCS=N elementary 
moves) and set J = 1. 

The typical behavior of the heat probability distributions (PD) is reported in Fig. [I] 
for the system with static coupling to the baths, T\ = 2.9, T 2 = 3, and a square geometry 
with L = M = 20. Much larger sizes can be hardly used because trajectories with a 
heat whose sign is opposite to that of the average value would be too rare. Results 
are qualitatively similar to those obtained with non-conserved dynamics [23]. Q^(r) 
(Q {t)) is on average negative (positive) and the relation 

(Q (1) (r)) + (Q( 2 )(r)) = (7) 

is verified. The PDs for Q} 1 ' and are similar but when AT = T 2 — T\ is increased 
the distribution of QP~> can be observed to be more peaked. 
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Figure 1. Heat PDs for (on the left) and (on the right) for a system evolving 
with spin-exchange-dynamics and T\ = 2.9, T 2 = 3, size 20x20 and t F = 5 • 10 s MCS. 
In the inset ai 1] V {Q {1) (r)) is plotted against {Q (1) (t) - (Q (1) (t))) / . Curves for 
different r collapse on a Gaussian mastercurve. The same collapse of data would be 
observed for Q( 2 \ 




Figure 2. Same parameters of Fig. [TJ log \V{Q, (n) {t)) /V{-Q (n) {t))] is plotted 
against (1/Ti — 1/T 2 )Q^(t) (n — 1 lower panel, n = 2 upper panel). 
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Figure 3. Same parameters of Fig. [TJ (different sizes). e(r, L) is plotted against t/L 
for different L. Results from a fit based on formula (fl4l) are also shown. 



The only characteristic time above T c is the (microscopic) relaxation time that can 
be measured from the decay of the autocorrelation function; it is of few MCS for the 
case of Fig. [TJ At sufficiently large r, greater than the relaxation time, due to the central 
limit theorem, one expects a Gaussian behavior for the probability distributions, namely 
V(Q^(r)) = ^(a^r^exph ^'^^ with {($>) ~ r and ~ ^F. 

This form is found with good accuracy, as shown by the collapse of the PD's at different 
times in the inset of Fig. [TJ 

In Fig. [2J in order to study the FR (j2j), the logarithm of the ratio 
^(Q (n) (r))/P(-QW(r)) is plotted as a function of A(3^ QW(r). For every value of 
t the data are well consistent with a linear relationship even if, for large values of the 
heat, the statistics becomes poor. The FR (jSJ) is verified if the slopes 

1ti 7>(Q<">(t)) 

n(»)/V\ = ^(-g ( "'(r)) /M 

1 J QW(r)A/3(") 1 J 

tend to 1 when r — > oo. In Fig. [3] the behavior of the distance e^ n ' = 1 — D^ n \r) from 
the expected asymptotic result is shown for the case of Fig. [TJ Indeed, this quantity 
goes to zero for large r showing the validity of Eq. (j2J). 

A similar behavior occurs in the case with dynamic couplings. In the upper panel 
of Fig. H]the PD's for the interface and bulk exchanged heats and are shown. 
The PD's corresponding to the colder bath are always higher and narrower. As in the 
case with static coupling, also now the PD's at different r can be rescaled on a single 
Gaussian master curve, as it can be seen in the inset of the upper panel of Fig. [H The 
distances tend to zero, as shown in Fig. [5j and the fluctuation relation (j2j) is verified. 
A scaling argument predicting the behavior of will be presented in Sec. 12.41 
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Figure 4. Upper panel (T 1 ,T 2 > T c ): Heat PDs for bulk heat £>W (on the left) 
and interface heat (on the right) for the system with dynamical coupling and 

Ti = 3.0, T 2 = 3.1, size 10x10 and t F = 5 ■ 10 8 MCS. In the inset ai^ V (Q w (r)) is 
plotted against (Q^ 1 ' (t) — (Q' 1 ' (t)))/ot • Curves for different r collapse on a Gaussian 
mastercurve. Lower panel (T%,T2 < T c ): Same kind of plot for T\ = 1,T 2 = 1.3, size 
10 x 10 and tp = 10 9 MCS. The inset shows that PD's at different r do not collapse 
on a single curve after rescaling. 



2.3. Deviations from the Green-Kubo relation 

As discussed in the introduction, once the validity of the FR has been ascertained in this 
context, it is interesting to see if deviations from the Green-Kubo relation are observed, 
in order to provide a rigid test of the FR not related to linear response theory. The GKR 
for the current J = \ < Q^(r) > |/r (which does not depend on n due to Eq. (JTj)) 
reads 

7 V 

l im 4^ = -T^> (9) 
at^o AT 2Tl v ; 
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where T> is related to the fluctuations of the heat Q exchanged with the bath in the 
equilibrium state at T\ through 

P= < e < T > 2 >. (10) 

T 

In order to check if the GKR ([9]) is verified we proceeded as follows: i) from an 
equilibrium simulation at T\ = 3 we have extracted T> and ii) by fixing T\ and varying 
AT in the range [0, 1] we have computed the current J . We have used the largest values 
of r for which the slope (jSJ) is measurable and the FR was verified. 

The data for the case with dynamical coupling to the baths are shown in Fig. [6j 
Our simulations clearly show deviations from the GKR ([9]) in a range of r where the 
FR holds. The same conclusion is arrived at in the case with static coupling to the 
reservoirs. This proves that in the cases considered insofar the FR is verified beyond 
the linear regime. 

2.4- Scaling behavior of finite-r corrections 

Generally, depend on the temperatures, the geometry of the system and on r. We 
propose now a scaling argument for the behavior of this quantity that takes also in 
account the nature of the coupling with the reservoirs. Considering a trajectory of 
length r in configuration space, enforcing Eq. ([6]), one can compute the ratio between 
the probability of the trajectory V(traj), conditioned to a given initial state, and that 
of the time-reversed evolution, obtaining 

e T i T 2 = e v ' T ™> , (11) 



V(-traj) 

where n' ^ n and AE = Qp^ (r) + (r) is the difference between the energies of the 
final and initial states. From Eq. after averaging over all possible trajectories 
it is straightforward to arrive to the following expression 



pstazt T \ 



e V) = ( e »¥»''g<«)M) i ( i2) 



where V staz (Q) (V staz (r)) is the probability of the initial (final) state. 

In order to proceed further, we make a quasi equilibrium hypothesis, namely 



j)stazf T \ 



e -/3iAE 1 -f3 2 AE 2 



with AE = AE\ + AE 2 . This is a strong assumption which is not expected to hold 
true in general. However it should be correct when the non-equilibrium drive is small, 
and the system is such that the NESS approaches the equilibrium state in the limit 
of small drive 0. In the present case, therefore, it is expected to apply for small AT, 

| The fulfillment of this condition is not obvious, particularly in systems with ergodicity breaking. It 
is not verified in the low temperature phase of the large-N model under shear flow [3T] , a model related 
(but with a mean field character) to the Ising one considered in this paper, nor in some driven mean 
field models, see 1521. 
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Figure 5. Same parameters of Fig. [4] (different sizes). Upper Panel (Ti,T2 > T c 
e^ n '(r, L) are plotted against r for different L. Lower panel (T\,T2 < T c ): : 
plotted against r for different L. 



which is indeed the case of our simulations, in the high temperature disordered phase. 
Furthermore, assuming that AE n and Q^ n ' are Gaussian distributed (which has been 
actually verified numerically for 7\,T2 > T c ), one finally arrives at 



\ 



(QW(r))\ 2 , v 2 



+ (A/JM)' (14) 



where the signs — , + are for n — 1,2, respectively, and v\ E is a quantity related to the 
variance of AE n which, for large r, is expected not to depend much on r. 

As discussed above, for large r, the quantities < Q^ n \r) > and (al^) 2 grow 
proportionally to r, so that their ratio is asymptotically constant. More precisely, 
since the heat is exchanged between neighboring spins coupled to different reservoirs, 
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indicating by Nfi ux the number of such couples of spins we expect Q 
Nfi ux T. In the model with static coupling to the baths one has Nfi ux oc L. With dynamic 
coupling, instead, since every spin in the system can feel one or the other temperature, 
one has Nfi ux oc N. Since v\ E is an extensive quantity proportional to the number N 
of spins, from Eq. (}Tlj) one obtains the scaling form e^(r, L) ~ f(x), with 



and f(x) is a (temperature dependent) scaling function with the large-x behavior 
f(x) ~ x^ 1 . Notice that the corrections to the universal law §Z§ are system dependent 
but with the same asymptotic dependence on r. We will now examine the validity of 
this scaling behavior in the cases considered so far. 

For the system with static coupling the data of Fig. [3] confirm our predictions: 
Curves with different L collapse when plotted against x = r/L, for all the values of x 
considered, and oc x~ l for sufficiently large r. We reported also the results obtained 
by applying formula ffT4l) with v\ E used as a fitting parameter; we observe a very good 
overlap with the numerical data. As expected, the fitting parameter results to be a 
quantity proportional to the total number N of lattice sites. Moreover, we considered 
different cases by varying Af3 in the interval 0.1-0.4 at fixed T%, and we did not find a 
significant dependence of the fitting parameter on Aj3. 

The forms ( fT4l) and ( fl5l) are also verified for the Ising model with dynamical coupling 
(see Fig. [5]). Sizes between L = 6 and L = 75 have been considered. The data collapse 
when plotted as a function of r (except, perhaps, for the system of size L = 6 which 
is probably too small for our scaling argument to fully apply). At large r one observes 
g( n ) r _1 , while preasymptotic results are again well reproduced by Eq. (Tl4"|) . 

2.5. Results for T 1; T 2 < T c 

Before discussing the results of the simulations in the low temperature phase it is 
useful to overview the behavior of the Ising model in contact with a single bath at 
temperature T < T c . In the thermodynamic limit the system is confined into one 
of the two pure states, which can be distinguished by the sign of the magnetization 
m(T) = (1/N) X^ili °i- This state is characterized by a finite relaxation time r eq (T). 
When N is finite, instead, genuine ergodicity breaking does not occur. Nevertheless, the 
system remains trapped into the basin of attraction of the pure states also in this case, 
although only for a finite time T erg (N,T). Therefore, for a finite size system in the low 
temperature phase, there is the additional timescale r erg (N, T), beside r eq (T), which can 
grow large. This ergodic time diverges when iV — > oo or T — ► 0, when r eq (T) is left as 
the only time scale in the system. For what follows, it is important to recall that, due 
to the trapping discussed above, an observation of a finite system on timescales much 
smaller than r erg (N,T) is representative of the state with ergodicity breaking which, 
however, strictly speaking, can be only realized for N = oo. 




static coupling to baths 
dynamic coupling to baths 



(15) 
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Figure 6. The quantity J / AT is plotted against AT for the system with dynamic 
coupling to the baths, size 10x10, T\ — 3. The horizontal dashed line is the value of 
V/(2Ti) obtained from the simulation of the system in the equilibrium state at T\. 
Near some of the points the values of the slopes (r) at r = 450 have been reported. 



A similar picture holds for systems in contact with two thermal baths, where there 
are two ergodic times r erg (N,T n ), one for each subsystem. Since the FR is expected to 
hold for t larger than the typical timescales of the system, it is interesting to study the 
role of the additional timescales r erg (N, T n ) on the FR and the interplay between r and 

We have studied numerically the model with dynamic coupling to the baths at 
T 1 ,T 2 < T c . T erg (N,T n ) can be evaluated as the time over which the autocorrelation 
functions C^{t - t') = (J2? =1 of \t)a^\t')), where crj n) denote spins in contact with 
the bath at T = T n , decay to zero (notice that r erg (N,Ti) > r erg (N, T 2 ) since Ti < T 2 ). 
Then, by varying T 1; T 2 and TV appropriately one can realize the limit of large r in the 
two cases with i) r r erg (N, T 2 ) or ii) r ^> r erg (N,Ti). In the former case, in the 
observation time-window r the system is practically confined into pure states while in 
the latter ergodicity is restored. Not surprisingly, in the latter case, we have observed 
a behavior very similar to that with Xi,T 2 > T c , except for the timescales required 
to access the asymptotic stage which are much larger in the low temperature regime, 
due to the tiny amount of heats exchanged. In particular, the PDs are found to be 
Gaussian. More interesting is the case i), for which the PDs are shown in the lower 
panel of Fig. HI The distributions extend on a smaller support, compared with those 
for Xi, T 2 > T c , even if the times over which the heats Q^ n ' are collected are larger (see 
insets of the upper panel of Fig. HJ). Moreover, quite interestingly, in this case the PDs 




Figure 7. The logarithm of the PDs of the lower panel of Fig. 0] (for only) 
(normalized by the height of A T of the peak) is plotted against \Q^(t) — (Q^ (r)\/a T . 
In the upper (lower) panel the branch QW < (qM) (q(i) > (Q«)) is plotted. The 
continuous line is the Gaussian behavior. In the inset the skewness is plotted against 
time. The dashed line is the power law behavior r -1 / 2 . 



are not Gaussian in the range of times accessed by the simulations. This can already 
be seen from the lower panel of Fig. HJ since the PDs are not symmetric with respect 
to the maximum. In particular, the right (left) tail of the distribution of (Q^) 
is visible fatter than the left (right) one. As a further test, in the inset of the same 
figure we tried to collapse the curves similarly to what done in the upper panel for the 
case above T c , but, due to the non Gaussian behavior, the collapse fails. In order to 
study the shape of the PDs more carefully, in Fig. [7] we plot the logarithm of the PD 
for (normalized by the height of A T of the peak) against \Q^ 1 \t) — (Q^(r))\/a T , 
where a\ is the variance of the distribution proportional to r. In this representation, 
a Gaussian PD should correspond to a power law, namely a straight line of slope 2. 
The figure shows that the Gaussian behavior is approached as r increases. For finite r, 
however, the distribution is strongly non Gaussian both around the average ~ (Q^ 1 ') 
and on the tails. In particular, the slope of the positive tail with ^> (Q^ 1 ') 
(lower panel) is always smaller than that of the negative one, confirming that the 
former is fatter. In addition, the approach to a Gaussian behavior is faster for the 
negative tail. This pattern of behavior can be expressed quantitatively by the skewness 
S(t) =< (Q(r)- < Q(t) >) 3 > / < (Q(t)- < Q(r) >) 2 > 3 / 2 of the distributions. 
This quantity is plotted in the inset of Fig. [3, showing a power law decay consistent 
with S(t) ~ r- a , with a ~ 1/2. 
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Regarding the slopes D^ n \r), they converge to 1 and the scaling e(r, L) ~ 1/x with 
x = r is verified, as shown in Fig. [5] (we have not reported the data for which are 
more noisy than those for Q^). We observe that a fit based on Eq. ( Ti~4"|) well reproduces 
the simulation data also in this case. As anticipated, the times required to reach the 
slope 1 are much longer than the corresponding ones at high temperature (but always 
smaller than r erg (N, T 2 ) > 5 x 10 8 MCS). These results suggest that the FR (T5]) and the 
scalings (I14|I15I) hold even in states with symmetry breaking and that the presence of the 
macroscopic timescales T erg (N,T n ) do not affect their validity. Recalling the discussion 
at the beginning of this section, it is reasonable to conjecture that for r <C r erg (N, T 2 ) 
the system enters a NESS which is representative of the state with broken symmetry that 
one would have (indefinitely) for an infinite system. For the same mechanism discussed 
above in equilibrium conditions, in this state there is a single relaxation time left (playing 
the role of T eq (T)), and the FR is obeyed on timescales larger then this microscopic time. 
Moreover, recalling the discussion of Sec. HJ this case is particularly interesting since the 
FR is verified (although with the preasymptotic corrections described by Eq. (1141) ) in a 
far from equilibrium situation in which the PDs are not Gaussian. 

3. Stationary states under shear 

In the present section we consider a system which is coupled to a single heat bath 
and is mechanically driven out of equilibrium. Mechanical work is done on the system 
by shearing it in the horizontal direction. Boundary conditions are periodic in the 
horizontal direction (parallel to the shear) and open in the vertical direction. In a shear 
event, the horizontal line with coordinate y G {1, • • • L} is moved by Xy lattice steps to 
the right. Shear events occur at regular intervals of r elementary Monte Carlo moves, 
with the shear period r submultiple of N = M x L. For each Monte Carlo sweep (MCS) 
over the lattice there are N/2 Kawasaki moves and N/r shear events. This model can 
be thought of as a discrete version for the kinetic equation of a binary mixture subject 
to the convective velocity v x (y) = Xy/{r/N) having a shear rate 7 = ^ = XN/r. The 
model, with single-spin-flip dynamics, was used to study domain growth properties of 
sheared systems in [33J. 

The heat released to the bath during the evolution will be computed as in Eq. ([5]), 
where the index (n) will be dropped since here the system is exchanging energy with a 
single reservoir at temperature T. Work done on the system is instead defined as the 
energy variation in shear events. 

Starting from a random configuration, with shear applied at a constant rate, a 
steady state is reached. As before, the PDs for heat and work relative to the steady 
state can be constructed collecting the values integrated over segments of length r in a 
long trajectory. 

In the following we report results obtained for a 50 x 4 lattice, with T = 5, shear 
step A = 1 and shear period r = 50, which corresponds to a shear rate 7 = 4. Data have 
been taken mainly with 10 7 MCS to reach the NESS and 10 8 MCS for the measurements. 
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100 200 



-Q(x), W(x) 

Figure 8. Heat and work PDs for the system under shear at different r. 5 • fO 7 
MCS have been discarded in order to reach stationarity and 5 ■ f 8 MCS have been 
used for measurements. See text for the other parameters. In the inset the curves are 
separately rescaled as in Fig. CD (<7 { t Q) and are the standard deviations of Q and 

W, respectively, for different t). 



A few longer (5 • 10 8 MCS) runs have been made for some data points. 

In Fig. [8] we show the work and heat PDs. At smaller r the heat distributions are 
more peaked than the work PDs while they become more similar when r increases. In 
all cases the PDs are very well fitted by Gaussian distributions and the work and heat 
PDs at different r collapse on a master curve when rescaled as in Fig. [IJ 

As in Sect. 2, also here we define the distance from the expected asymptotic 
behavior e^°> = 1 — D^°\r) where, analogously to Eq. ([8]), the slope D^°\r) is defined 
as 

d(0){t) = (16) 

where the signs +, — are for O = W, Q, respectively. In Fig. [9] we report the behavior 
of as a function of r. tend to zero, but differently than in the case of 

thermodynamic driving, the approach to the expected asymptotic result is slower, 
especially for the heat. 

Analogously to the case of two thermostats, the ratio of the probability of a 
trajectory over that of the time- reversed is given by 

Vitraj) Q(r) W(t) AE 

,-— = e ~ 5-, (17) 



V(-traj) 

where, for each trajectory, the relation Q(t) + W(r) = AE holds. Proceeding as in the 
case with thermodynamic drive one arrives at an expression similar to Eq. ([1] 



(g) _ T(Q(r)) 



\ 



.Jo ))2 + {a io Y W 



where the signs — , + are for W, Q respectively, and v 2 is an unknown quantity that 
can be used as a fit parameter. In doing that we find a very good agreement between 
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Figure 9. e' ' (O = Q, W) as a function of r for the same system of Fig. [8] 



formula (|T8l) and the numerical data for IV, while the accordance (not shown in Fig. [9]) 
is very poor for Q, probably because, as already noticed, this quantity has a slower 
convergence. The faster convergence of W can be possibly ascribed to the fact that, 
after making the quasi-equilibrium assumption Vstaz ^ — e~ p and calculating the 
averages over trajectories as in Sec. 12.41 the boundary term proportional to AE cancels 
for the work while it remains when one consider the behavior of the heat. 

4. Discussion and Conclusions 

In this Article we have considered the issue of the FRs in the Ising model with nearest 
neighbor interactions. We focused on two different NESS, induced by a thermodynamic 
drive provided by two thermostats at different temperatures, or by a mechanical forcing 
realized by shearing the system. We find that the FR (QQ) is verified in any case when 
t is sufficiently large. For the case with different thermostats, we have also provided 
an exacting test of the FR in the far from equilibrium regime where deviations from 
the linear response theory are observed with the Green-Kubo relation violated (above 
T c ), or the PDs are not Gaussian (below T c ). Furthermore, we have analyzed the effects 
of finite time by proposing a scaling law which takes into account the geometry of 
the system and the details of the interaction with the baths and that is well verified 
numerically. According to this law the leading corrections to the FR ([1]) decay as r -1 , 
in agreement with what found in some previous studies [3], [TT], [TTj , while the dependence 
on the system size relies on the details of the coupling to the baths. The derivation of 
our scaling law is based on few strong but general hypotheses, which could, in principle, 
be verified also in different systems. A natural perspective for future work, therefore, is 
the understanding of the possible generality of our results in different contexts. 

Most of the results of this paper are obtained for systems that are able to equilibrate 
in the small entropy production limit [32]. Basically this means that, for this class of 
systems, the equilibrium state is recovered in a finite time when the drive is switched off. 
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In this case the FR generally holds, as we have explicitly verified. The generalization 
of the FR to systems which do not equilibrate in the small entropy production limit 
is a very interesting issue. Non-equilibration typically happens when the limit of large 
t is taken after the thermodynamic limit in systems, such as glasses or ferromagnets 
below the critical temperature, whose equilibration time diverges with the system size. 
Regarding these non-equilibrating systems, some conjectures have recently appeared 
in the literature [31], inferring that a relation similar to (Tj[|) holds, where, however, 
the entropy production, instead of being a function of the bath(s) temperature (as for 
instance in Eqs. (l2|3|) ). depends on T e ff, namely the effective temperature [35] which 
characterizes the dynamics of these systems. For the model considered in this paper one 
has T e ff = oo [36]. In this perspective, therefore, it would be very interesting to study 
the technically demanding case of the Ising model with shear for T < T c . 
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